N<-101	#入力次元
JRandSd<- sqrt(1)	#標準偏差=sqrt(分散)
JRandM<-0	#平均
XRandSd<- 1/sqrt(N)	#標準偏差
XRandM<-0	#平均
eta<- .5
M<- N/eta
MML<- 2.5*N/eta

#set.seed(300)
b<-matrix(rnorm(N,JRandM,JRandSd),N,1)
b<-b*sqrt(N/sum(b^2))
bL<-sqrt(sum(b^2))
#set.seed(500)
j1<-matrix(rnorm(N,JRandM,JRandSd),N,1)
j1<-j1*sqrt(N/sum(j1^2))
#set.seed(800)
j2<-matrix(rnorm(N,JRandM,JRandSd),N,1)
j2<-j2*sqrt(N/sum(j2^2))

R<-matrix(0,MML+M+1,2)
q<-matrix(0,MML+M+1,1)
L<-matrix(0,MML+M+1,2)

tlearn <- function(){
	for(i in 1:M) {
		if(i==1||((R[(i-1),1]) < 0.6)){
			x<-matrix(rnorm(N,mean=XRandM,sd=XRandSd),N,1)
			ub<-sum(x*b)
			uj1<-sum(x*j1)
			j1 <<- j1 + eta*(ub-uj1)*x
		}
		
		if(i==1||((R[(i-1),2]) < 0.6)){
			x<-matrix(rnorm(N,mean=XRandM,sd=XRandSd),N,1)
			ub<-sum(x*b)
			uj2<-sum(x*j2)
			j2 <<- j2 + eta*(ub-uj2)*x
		}
		
		L[i,1]<<-sqrt(sum(j1^2))
		R[i,1]<<-sum(j1*b)/(L[i,1]*bL)
		L[i,2]<<-sqrt(sum(j2^2))
		R[i,2]<<-sum(j2*b)/(L[i,2]*bL)
		q[i]<<-sum(j1*j2)/(L[i,1]*L[i,2])
		
		if(i!=1&&(0.6<(R[(i-1),1]))&&(0.6<(R[(i-1),2]))){
			M<<-i
			break()
		}
	}
}
mlearn <- function(){
	for (i in (M+1):(M+MML+1)) {
		x<-matrix(rnorm(N,XRandM,XRandSd),N,1)
		uj1<-sum(x*j1)
		uj2<-sum(x*j2)
		j1 <<- j1 + eta*(uj2-uj1)*x
		j2 <<- j2 + eta*(uj1-uj2)*x
		L[i,1]<<-sqrt(sum(j1^2))
		L[i,2]<<-sqrt(sum(j2^2))
		R[i,1]<<-sum(j1*b)/(L[i,1]*bL)
		R[i,2]<<-sum(j2*b)/(L[i,2]*bL)
		q[i]<<-sum(j1*j2)/(L[i,1]*L[i,2])
	}
}
#MAIN routine
tlearn()
j1<-j1*sqrt(N/sum(j1^2))
j2<-j2*sqrt(N/sum(j2^2))
mlearn()
#
plot(L[(M+1):(M+MML+1),1]/sqrt(N),R[(M+1):(M+MML+1),1],xlim=c(0.75,1),ylim=c(0.5,0.8),type="l")
par(new=T);plot(L[(M+1):(M+MML+1),2]/sqrt(N),R[(M+1):(M+MML+1),2],xlim=c(0.75,1),ylim=c(0.5,0.8),type="l")
#
plot(L[1:(M+MML+1),1]/sqrt(N),R[1:(M+MML+1),1],type="l")
par(new=T);plot(L[1:(M+MML+1),2]/sqrt(N),R[1:(M+MML+1),2],type="l")
#
plot(1:(MML+M+1),R[1:(MML+M+1),1],ylim=c(0,1),type="l")
par(new=T);plot(1:(MML+M+1),R[1:(MML+M+1),2],ylim=c(0,1),type="l")
par(new=T);plot(1:(MML+M+1),q[1:(MML+M+1)],ylim=c(0,1),type="l")
#
plot(M:(MML+M+1),R[M:(MML+M+1),1],ylim=c(0,1),type="l")
par(new=T);plot(M:(MML+M+1),R[M:(MML+M+1),2],ylim=c(0,1),type="l")
par(new=T);plot(M:(MML+M+1),q[M:(MML+M+1)],ylim=c(0,1),type="l")
